This is R Markdown file produced to evidence the spatial data analysis skills learned throughout CASA0005: Geographical Information System and Sience, specifically for part 2 of its assessment.
setwd("D:/Xin/MSc/Modules/CASA0005_GISS/Assessment/Assessment Part 2/")
library(geojsonio)
##
## Attaching package: 'geojsonio'
## The following object is masked from 'package:base':
##
## pretty
library(tmap)
hunt <- geojson_read("https://www.dropbox.com/s/wa2ip35tcmt93g3/Team7.geojson?raw=1", what = "sp")
tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(hunt) + tm_lines(col="red", lwd=4)
## Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
##########################################
####Read in tube stations
library(sf)
stn <- st_read("https://www.doogal.co.uk/LondonStationsKML.ashx")
## Reading layer `London stations with zone information' from data source `https://www.doogal.co.uk/LondonStationsKML.ashx' using driver `KML'
## Simple feature collection with 641 features and 2 fields
## geometry type: POINT
## dimension: XYZ
## bbox: xmin: -0.61143 ymin: 51.28216 xmax: 0.329952 ymax: 51.74702
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defs
ukgrid="+init=epsg:27700"
stn<-st_transform(stn,ukgrid)
tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(stn)+tm_markers(col="red",size=0.01,shape=10)
## Symbol shapes other than circles or icons are not supported in view mode.
##########################################
####geocode treasure locations
library(ggmap)
## Loading required package: ggplot2
## Google Maps API Terms of Service: http://developers.google.com/maps/terms.
## Please cite ggmap if you use it: see citation("ggmap") for details.
library(tidyverse)
## -- Attaching packages ---------------------------------------------------------------------------------- tidyverse 1.2.1 --
## v tibble 1.4.2 v purrr 0.2.5
## v tidyr 0.8.1 v dplyr 0.7.7
## v readr 1.1.1 v stringr 1.3.1
## v tibble 1.4.2 v forcats 0.3.0
## -- Conflicts ------------------------------------------------------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
library(RCurl)
## Loading required package: bitops
##
## Attaching package: 'RCurl'
## The following object is masked from 'package:tidyr':
##
## complete
library(RJSONIO)
## Warning: package 'RJSONIO' was built under R version 3.5.2
library(plyr)
## -------------------------------------------------------------------------
## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)
## -------------------------------------------------------------------------
##
## Attaching package: 'plyr'
## The following objects are masked from 'package:dplyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarize
## The following object is masked from 'package:purrr':
##
## compact
#Read in the list of treasure locations and their points scored from the dropbox given
huntaddresses <- read_csv("https://www.dropbox.com/s/v66l4cx7aia9jlo/Locations.csv?raw=1")
## Parsed with column specification:
## cols(
## Location = col_character(),
## Points = col_integer()
## )
#Set google map api for geocoding purposes
my_api="AIzaSyA5EUX3ZcX4u6b_OhogTM9FaWlxu4U5fJc"
register_google(key=my_api)
ggmap(get_googlemap())
## Source : https://maps.googleapis.com/maps/api/staticmap?center=29.763284,-95.363271&zoom=10&size=640x640&scale=2&maptype=terrain&key=AIzaSyA5EUX3ZcX4u6b_OhogTM9FaWlxu4U5fJc
#highlight this whole block and create this function to access the Google Places API
url <- function(address, return.call = "json", sensor = "false") {
root <- "https://maps.google.com/maps/api/geocode/"
u <- paste(root, return.call, "?address=", address, "&sensor=", sensor, "&key=", my_api, sep = "")
return(URLencode(u))
}
#highlight this whole block and create this function to geocode some places just from a random list of treasure hunt locations
geoCode <- function(address,verbose=FALSE) {
if(verbose) cat(address,"\n")
u <- url(address)
doc <- getURL(u)
x <- fromJSON(doc,simplify = FALSE)
if(x$status=="OK") {
lat <- x$results[[1]]$geometry$location$lat
lng <- x$results[[1]]$geometry$location$lng
location_type <- x$results[[1]]$geometry$location_type
formatted_address <- x$results[[1]]$formatted_address
return(c(lat, lng, location_type, formatted_address))
Sys.sleep(0.5)
} else {
return(c(NA,NA,NA, NA))
}
}
#now use the geoCode() function (which calls the URL function) to geocode our list of places
#for loop to cycle through every treasure hunt location
i=1
for(i in 1:nrow(huntaddresses)){
# Every nine records, pause 3 seconds so that the API doesn't kick us off...
if(i %% 9 == 0) Sys.sleep(3)
#now create a temporary list of useful elements
tempdf <- as.list(geoCode(huntaddresses[i,1]))
#and write these back into our dataframe
huntaddresses[i,3] <- tempdf[1]
huntaddresses[i,4] <- tempdf[2]
huntaddresses[i,5] <- tempdf[4]
}
# rename the columns
names(huntaddresses) <- c("Location","Points","lat","lon","GoogleAddress")
head(huntaddresses)
## # A tibble: 6 x 5
## Location Points lat lon GoogleAddress
## <chr> <int> <chr> <chr> <chr>
## 1 Picadilly Circus, ~ 2 51.509~ -0.13~ Piccadilly Circus, London W1B~
## 2 Nelson's Column, L~ 2 51.507~ -0.12~ 5 Trafalgar Square, London WC~
## 3 Big Ben, London 2 51.500~ -0.12~ Westminster, London SW1A 0AA,~
## 4 100 Club, Oxford S~ 5 51.516~ -0.13~ Century House, 100 Oxford St,~
## 5 Fabric, Charterhou~ 5 51.519~ -0.10~ 77A Charterhouse St, Clerkenw~
## 6 The Gherkin, 30 St~ 2 51.514~ -0.08~ 30 St Mary Axe, London EC3A 8~
#write a new .csv file to your working directory
write_csv(huntaddresses, "huntaddresses.csv")
##########################################
####Transform hunt data into simple feature objects, and measure the distance travelled on the hunt route
hunt_sf<-st_as_sf(hunt,ukgrid)
st_length(hunt_sf)
## 46610.68 [m]
#Transform hunt data into simple feature objects
hunt_sf<-st_transform(hunt_sf,ukgrid)
#Plot hunt route and stations on the same map
tm_shape(stn)+tm_markers(col="red",size=0.02,shape=16)+tm_shape(hunt_sf)+tm_lines(col="navy",lwd=4)
##########################################
####Create station buffers of 100m radii
stn_buff<-st_buffer(stn,dist=100)
#View station buffers and hunt route on the same map
tm_shape(stn_buff)+tm_polygons(col="red",size=0.02,shape=16,alpha=0.3)+tm_shape(hunt_sf)+tm_lines(col="navy",lwd=4)
#Count the instances of intersects between the station buffers created and the hunt route.
num_stn_intersects<-sum(lengths(st_intersects(stn_buff,hunt_sf)))
num_stn_intersects
## [1] 24
#Set latitude and longitude coordinate system
latlon<-"+init=epsg:4326"
#Retrieve coordinates data from huntaddresses data generated previously
huntaddresses$lat<-as.numeric(as.character(huntaddresses$lat))
huntaddresses$lon<-as.numeric(as.character(huntaddresses$lon))
#Create simple features from huntaddresses dataframe
huntadd_sf<-st_as_sf(x=huntaddresses,coords=c("lon","lat"),crs=latlon)
#And convert fro 4326 to 27700 crs
huntadd_sf<-st_transform(huntadd_sf,ukgrid)
#Check crs
st_crs(huntadd_sf)
## Coordinate Reference System:
## EPSG: 27700
## proj4string: "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=446.448,-125.157,542.06,0.1502,0.247,0.8421,-20.4894 +units=m +no_defs"
#Create buffers for huntaddresses
huntadd_buff<-st_buffer(huntadd_sf,dist=300)
#View hunt location buffers and hunt route on one map
tm_shape(huntadd_buff)+tm_polygons(col="red",size=0.02,shape=16,alpha=0.3)+tm_shape(hunt_sf)+tm_lines(col="navy",lwd=4)
#Count the instances of intersects between the huntaddress buffers and the hunt route
num_huntadd_intersects<-sum(lengths(st_intersects(huntadd_buff,hunt_sf)))
num_huntadd_intersects
## [1] 16
##########################################
####Calculate max and min male life expectancy
#Read in the England shapefile as sp dataframe
library(geojsonio)
EnglandMap<-geojson_read("http://geoportal.statistics.gov.uk/datasets/8edafbe3276d4b56aec60991cbddda50_2.geojson", what = "sp")
#Crop the London shapefile only
LondonMap <- EnglandMap[grep("^E09",EnglandMap@data$lad15cd),]
#Quickly plot the London wards shapefile
qtm(LondonMap)
#Convert sp objects into simple features
LondonMap_sf<-st_as_sf(LondonMap)
#Convert the crs to 27700 and check
LondonMap_sf<-st_transform(LondonMap_sf,ukgrid)
st_crs(LondonMap_sf)
## Coordinate Reference System:
## EPSG: 27700
## proj4string: "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=446.448,-125.157,542.06,0.1502,0.247,0.8421,-20.4894 +units=m +no_defs"
#Read in male life expectancy data
library(readxl)
mle_raw<-data.frame(read_excel("leatbirthandatage65byukla201517.xls",sheet=2,col_names=TRUE,col_types=NULL,na="",skip=2))
library(tmaptools)
LondonMap_sf<-append_data(LondonMap_sf,mle_raw,key.shp="lad15cd",key.data="Area.Codes",ignore.duplicates=TRUE,ignore.na=TRUE)
## Data contains duplicated keys: NA
## mle_raw key variable "Area.Codes" contains NA's, which are ignored
## Over coverage: 411 out of 444 data records were not appended. Run over_coverage() to get the corresponding data row numbers and key values.
#Append indicator column for intersect
LondonMap_sf$IntersectIndicator<-0
#Assign 1 as indicator where wards intersect with hunt route
for(i in 1:nrow(LondonMap_sf)) {
if(sum(lengths(st_intersects(LondonMap_sf[i,],hunt_sf)))==1) {
LondonMap_sf[i,ncol(LondonMap_sf)]=1
} else {
LondonMap_sf[i,ncol(LondonMap_sf)]=0
}
}
#Slice the London wards where hunt route intersect
subset<-subset(LondonMap_sf,LondonMap_sf$IntersectIndicator==1)
#The intersecting London ward with the highest male life expectancy in 2015-2017 at births, and the years of expectancy
subset[which.min(subset$X2015.2017),]$lad15nm
## [1] Lambeth
## 380 Levels: Aberdeen City Aberdeenshire Adur Allerdale ... York
max(subset$X2015.2017,na.rm=T)
## [1] 82.66597
#The intersecting London ward with the lowest male life expectancy in 2015-2017 at births, and the years of expectancy
subset[which.max(subset$X2015.2017),]$lad15nm
## [1] Westminster
## 380 Levels: Aberdeen City Aberdeenshire Adur Allerdale ... York
min(subset$X2015.2017,na.rm=T)
## [1] 78.69568
#Calculate the average years of male life expectancies in 2015-2017 at births
sum(subset$X2015.2017,na.rm = T)/nrow(subset)
## [1] 71.08626
##########################################
####Test spatial randomness of hunt locations
#create London boundaries with Northing/Easting, and set a window as the London boundary
library(spatstat)
## Warning: package 'spatstat' was built under R version 3.5.2
## Loading required package: spatstat.data
## Warning: package 'spatstat.data' was built under R version 3.5.2
## Loading required package: nlme
##
## Attaching package: 'nlme'
## The following object is masked from 'package:dplyr':
##
## collapse
## Loading required package: rpart
##
## spatstat 1.57-1 (nickname: 'Cartoon Physics')
## For an introduction to spatstat, type 'beginner'
library(rgdal)
## Warning: package 'rgdal' was built under R version 3.5.2
## Loading required package: sp
## rgdal: version: 1.3-6, (SVN revision 773)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
## Path to GDAL shared files: C:/Users/Xin/Documents/R/win-library/3.5/sf/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
## Path to PROJ.4 shared files: C:/Users/Xin/Documents/R/win-library/3.5/sf/proj
## Linking to sp version: 1.3-1
LondonMapBNG<-spTransform(LondonMap,ukgrid)
summary(LondonMapBNG)
## Object of class SpatialPolygonsDataFrame
## Coordinates:
## min max
## x 503576.3 561958.7
## y 155850.7 200934.0
## Is projected: TRUE
## proj4string :
## [+init=epsg:27700 +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717
## +x_0=400000 +y_0=-100000 +datum=OSGB36 +units=m +no_defs
## +ellps=airy
## +towgs84=446.448,-125.157,542.060,0.1502,0.2470,0.8421,-20.4894]
## Data attributes:
## lad15cd lad15nm lad15nmw objectid
## E09000001: 1 Barking and Dagenham: 1 :33 Min. :294
## E09000002: 1 Barnet : 1 Abertawe : 0 1st Qu.:302
## E09000003: 1 Bexley : 1 Blaenau Gwent: 0 Median :310
## E09000004: 1 Brent : 1 Bro Morgannwg: 0 Mean :310
## E09000005: 1 Bromley : 1 Caerdydd : 0 3rd Qu.:318
## E09000006: 1 Camden : 1 Caerffili : 0 Max. :326
## (Other) :27 (Other) :27 (Other) : 0
## st_lengthshape st_areashape
## Min. : 8929 Min. : 2897649
## 1st Qu.:28384 1st Qu.: 26797942
## Median :37664 Median : 37628571
## Mean :39255 Mean : 47682317
## 3rd Qu.:46679 3rd Qu.: 56413925
## Max. :74641 Max. :150125298
##
window<-as.owin(LondonMapBNG)
#Create a point pattern object for the hunt locations
huntadd.ppp<-ppp(x=st_coordinates(huntadd_sf)[,1],y=st_coordinates(huntadd_sf)[,2],window=window)
## Warning in ppp(x = st_coordinates(huntadd_sf)[, 1], y =
## st_coordinates(huntadd_sf)[, : 1 out of 50 points had NA or NaN coordinate
## values, and was discarded
## Warning: 2 points were rejected as lying outside the specified window
huntadd.ppp<-as.ppp(huntadd.ppp)
plot(huntadd.ppp,pch=16,cex=1,main="Hunt Locations in London")
#Plot the kernel density estimation map
plot(density(huntadd.ppp, sigma = 1000))
#Run a simple Ripley's K test on the hunt location data in London boundar, the maximum radius plotted is set at 5km.
K<-Kest(huntadd.ppp,correction="border",rmax=5000)
plot(K)